global n k casenum kmax draws
% This file conducts the GMM IV-Logit estimation on turnout for the PRI,
% and the opposition taking abstention as baseline
% Tobias Pfutze, May 2009

% Upload Data in order: pripart oppmaxpart remit_prophh_intl
% dist_border cycleB1 cycleB2 zona_norte labinc_perhh_log
% labinc_perhh_stdev illiterate indig margmuni00 poptot_2000 unempl
% raildistjuarez_05 histmig24_edoprop00 oppwin_since80 pripart1 oppmaxpart1

RAW=load('/Users/tpfutze/Dropbox/Research/Participation/Matlab/Data/participation_adjall.txt');
[N,M]=size(RAW);

% Replace total population by its logged value
RAW(:,13)=log(RAW(:,13));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Start by creating single random vector for bootstrap %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
draws=10000;
rand('state',1812);
unifrv=rand(draws,1);

% Create two separate matrices conditional on whether municipality has 
% had an opposition party in office since 1980
s0=1;
s1=1;
for i=1:N
    if (RAW(i,17)==0)
        RAW0(s0,1:M)=RAW(i,1:M);
        s0=s0+1;
    elseif RAW(i,17)==1
        RAW1(s1,1:M)=RAW(i,1:M);
        s1=s1+1;
    end
end
[N0,M0]=size(RAW0);
[N1,M1]=size(RAW1);


% Create vectors with randomly selected observations for bootstrap
RAND=ceil(N*unifrv);
RAND0=ceil(N0*unifrv);
RAND1=ceil(N1*unifrv);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create all matrices of dep., indep. and instr. variables %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create matrices and vectors 

%Dependent variable are log odds ratios
y=log(RAW(:,1:2)./(kron(ones(1,2),(1-RAW(:,1)-RAW(:,2)))));
y0=log(RAW0(:,1:2)./(kron(ones(1,2),(1-RAW0(:,1)-RAW0(:,2)))));       
y1=log(RAW1(:,1:2)./(kron(ones(1,2),(1-RAW1(:,1)-RAW1(:,2)))));

x1a=[ones(N,1) RAW(:,3:7) RAW(:,18)];                               
x2a=[ones(N,1) RAW(:,3:14) RAW(:,18)];
x1b=[ones(N,1) RAW(:,3:7) RAW(:,19)];                               
x2b=[ones(N,1) RAW(:,3:14) RAW(:,19)];

x01a=[ones(N0,1) RAW0(:,3:7) RAW0(:,18)];                               
x02a=[ones(N0,1) RAW0(:,3:14) RAW0(:,18)]; 
x01b=[ones(N0,1) RAW0(:,3:7) RAW0(:,19)];                               
x02b=[ones(N0,1) RAW0(:,3:14) RAW0(:,19)];

x11a=[ones(N1,1) RAW1(:,3:7) RAW1(:,18)];                               
x12a=[ones(N1,1) RAW1(:,3:14) RAW1(:,18)]; 
x11b=[ones(N1,1) RAW1(:,3:7) RAW1(:,19)];                               
x12b=[ones(N1,1) RAW1(:,3:14) RAW1(:,19)]; 

z1a=[ones(N,1) RAW(:,15:16) RAW(:,4:7) RAW(:,18)];                               
z2a=[ones(N,1) RAW(:,15:16) RAW(:,4:14) RAW(:,18)];
z1b=[ones(N,1) RAW(:,15:16) RAW(:,4:7) RAW(:,19)];                               
z2b=[ones(N,1) RAW(:,15:16) RAW(:,4:14) RAW(:,19)];

z01a=[ones(N0,1) RAW0(:,15:16) RAW0(:,4:7) RAW0(:,18)];                               
z02a=[ones(N0,1) RAW0(:,15:16) RAW0(:,4:14) RAW0(:,18)];
z01b=[ones(N0,1) RAW0(:,15:16) RAW0(:,4:7) RAW0(:,19)];                               
z02b=[ones(N0,1) RAW0(:,15:16) RAW0(:,4:14) RAW0(:,19)];

z11a=[ones(N1,1) RAW1(:,15:16) RAW1(:,4:7) RAW1(:,18)];                               
z12a=[ones(N1,1) RAW1(:,15:16) RAW1(:,4:14) RAW1(:,18)]; 
z11b=[ones(N1,1) RAW1(:,15:16) RAW1(:,4:7) RAW1(:,19)];                               
z12b=[ones(N1,1) RAW1(:,15:16) RAW1(:,4:14) RAW1(:,19)];



%%%%%%%%%%%%%%%%%%%
% DEFINE MATRICES %
%%%%%%%%%%%%%%%%%%%

MEAN1a=[mean(x1a(:,1:3)) 0 0 0 mean(x1a(:,7))];
MEAN1b=[mean(x1b(:,1:3)) 0 0 0 mean(x1b(:,7))];
MEAN2a=[mean(x2a(:,1:3)) 0 0 0  mean(x2a(:,7:14))];
MEAN2b=[mean(x2b(:,1:3)) 0 0 0  mean(x2b(:,7:14))];


%Matrices for point estimates and significance tests
[nmax,kmax]=size(x2a);

BETA=zeros(2*kmax,6);    % betas for shares
F_WALD=zeros(3,6);
T_WALD=zeros(2*kmax,6);
TJ_WALD=zeros(kmax,6);
PARTIAL=zeros(2*(kmax-1),6);
DOF=zeros(2,6);
OIR=zeros(1,6);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%% PARTICIPATION RESULTS %%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%% ENTIRE SAMPLE %%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1ST CASE
casenum=1;
[n,k] = size(x1a);
[n,l] = size(z1a);
mean_a = MEAN1a;
mean_b = MEAN1b;
Y=(reshape(y',1,2*n))';
for i=1:n
    X((i-1)*2+1:(i-1)*2+2,1:2*k)=[kron([1;0],x1a(i,:)) kron([0;1],x1b(i,:))];
    Z((i-1)*2+1:(i-1)*2+2,1:2*l)=[kron([1;0],z1a(i,:)) kron([0;1],z1b(i,:))];
end


[beta, T, T_J, F, part, oir] = SURPartiesRobustLagvar_IVGMM(Y, X, Z, mean_a, mean_b, RAND, n, k, l);


% Put into final matrices
BETA(1:k, casenum)=beta(1:k);
BETA(kmax+1:kmax+k, casenum)=beta(k+1:2*k);
T_WALD(1:k, casenum)=T(1:k);
T_WALD(kmax+1:kmax+k, casenum)=T(k+1:2*k);
TJ_WALD(1:k, casenum)=T_J(1:k);
F_WALD(:, casenum)=F;
PARTIAL(1:k-1, casenum)=part(1:k-1);
PARTIAL(kmax:kmax+(k-2), casenum)=part(k:2*k-2);
DOF(:,casenum)=[k, l];
OIR(:, casenum)=oir;
clear Y X Z P beta  W F T T_J alpha ETA part mean mean_l mean_h oir;
 


% 2ND CASE
casenum=casenum+1;
[n,k] = size(x2a);
[n,l] = size(z2a);
mean_a = MEAN2a;
mean_b = MEAN2b;
Y=(reshape(y',1,2*n))';
for i=1:n
    X((i-1)*2+1:(i-1)*2+2,1:2*k)=[kron([1;0],x2a(i,:)) kron([0;1],x2b(i,:))];
    Z((i-1)*2+1:(i-1)*2+2,1:2*l)=[kron([1;0],z2a(i,:)) kron([0;1],z2b(i,:))];
end


[beta, T, T_J, F, part, oir] = SURPartiesRobustLagvar_IVGMM(Y, X, Z, mean_a, mean_b, RAND, n, k, l);


% Put into final matrices
BETA(1:k, casenum)=beta(1:k);
BETA(kmax+1:kmax+k, casenum)=beta(k+1:2*k);
T_WALD(1:k, casenum)=T(1:k);
T_WALD(kmax+1:kmax+k, casenum)=T(k+1:2*k);
TJ_WALD(1:k, casenum)=T_J(1:k);
F_WALD(:, casenum)=F;
PARTIAL(1:k-1, casenum)=part(1:k-1);
PARTIAL(kmax:kmax+(k-2), casenum)=part(k:2*k-2);
DOF(:,casenum)=[k, l];
OIR(:, casenum)=oir;
clear Y X Z P beta  W F T T_J alpha ETA part mean mean_l mean_h oir;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%   ALWAYS PRI   %%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1ST CASE
casenum=casenum+1;
[n,k] = size(x01a);
[n,l] = size(z01a);
mean_a = MEAN1a;
mean_b = MEAN1b;
Y=(reshape(y0',1,2*n))';
for i=1:n
    X((i-1)*2+1:(i-1)*2+2,1:2*k)=[kron([1;0],x01a(i,:)) kron([0;1],x01b(i,:))];
    Z((i-1)*2+1:(i-1)*2+2,1:2*l)=[kron([1;0],z01a(i,:)) kron([0;1],z01b(i,:))];
end

[beta, T, T_J, F, part, oir] = SURPartiesRobustLagvar_IVGMM(Y, X, Z, mean_a, mean_b, RAND0, n, k, l);


% Put into final matrices
BETA(1:k, casenum)=beta(1:k);
BETA(kmax+1:kmax+k, casenum)=beta(k+1:2*k);
T_WALD(1:k, casenum)=T(1:k);
T_WALD(kmax+1:kmax+k, casenum)=T(k+1:2*k);
TJ_WALD(1:k, casenum)=T_J(1:k);
F_WALD(:, casenum)=F;
PARTIAL(1:k-1, casenum)=part(1:k-1);
PARTIAL(kmax:kmax+(k-2), casenum)=part(k:2*k-2);
DOF(:,casenum)=[k, l];
OIR(:, casenum)=oir;
clear Y X Z P beta  W F T T_J alpha ETA part mean mean_l mean_h oir;
 


% 2ND CASE
casenum=casenum+1;
[n,k] = size(x02a);
[n,l] = size(z02a);
mean_a = MEAN2a;
mean_b = MEAN2b;
Y=(reshape(y0',1,2*n))';
for i=1:n
    X((i-1)*2+1:(i-1)*2+2,1:2*k)=[kron([1;0],x02a(i,:)) kron([0;1],x02b(i,:))];
    Z((i-1)*2+1:(i-1)*2+2,1:2*l)=[kron([1;0],z02a(i,:)) kron([0;1],z02b(i,:))];
end


[beta, T, T_J, F, part, oir] = SURPartiesRobustLagvar_IVGMM(Y, X, Z, mean_a, mean_b, RAND0, n, k, l);


% Put into final matrices
BETA(1:k, casenum)=beta(1:k);
BETA(kmax+1:kmax+k, casenum)=beta(k+1:2*k);
T_WALD(1:k, casenum)=T(1:k);
T_WALD(kmax+1:kmax+k, casenum)=T(k+1:2*k);
TJ_WALD(1:k, casenum)=T_J(1:k);
F_WALD(:, casenum)=F;
PARTIAL(1:k-1, casenum)=part(1:k-1);
PARTIAL(kmax:kmax+(k-2), casenum)=part(k:2*k-2);
DOF(:,casenum)=[k, l];
OIR(:, casenum)=oir;
clear Y X Z P beta  W F T T_J alpha ETA part mean mean_l mean_h oir;




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%% ALREADY OPPOSITION %%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1ST CASE
casenum=casenum+1;
[n,k] = size(x11a);
[n,l] = size(z11a);
mean_a = MEAN1a;
mean_b = MEAN1b;
Y=(reshape(y1',1,2*n))';
for i=1:n
    X((i-1)*2+1:(i-1)*2+2,1:2*k)=[kron([1;0],x11a(i,:)) kron([0;1],x11b(i,:))];
    Z((i-1)*2+1:(i-1)*2+2,1:2*l)=[kron([1;0],z11a(i,:)) kron([0;1],z11b(i,:))];
end


[beta, T, T_J, F, part, oir] = SURPartiesRobustLagvar_IVGMM(Y, X, Z, mean_a, mean_b, RAND1, n, k, l);


% Put into final matrices
BETA(1:k, casenum)=beta(1:k);
BETA(kmax+1:kmax+k, casenum)=beta(k+1:2*k);
T_WALD(1:k, casenum)=T(1:k);
T_WALD(kmax+1:kmax+k, casenum)=T(k+1:2*k);
TJ_WALD(1:k, casenum)=T_J(1:k);
F_WALD(:, casenum)=F;
PARTIAL(1:k-1, casenum)=part(1:k-1);
PARTIAL(kmax:kmax+(k-2), casenum)=part(k:2*k-2);
DOF(:,casenum)=[k, l];
OIR(:, casenum)=oir;
clear Y X Z P beta  W F T T_J alpha ETA part mean mean_l mean_h oir;
 


% 2ND CASE
casenum=casenum+1;
[n,k] = size(x12a);
[n,l] = size(z12a);
mean_a = MEAN2a;
mean_b = MEAN2b;
Y=(reshape(y1',1,2*n))';
for i=1:n
    X((i-1)*2+1:(i-1)*2+2,1:2*k)=[kron([1;0],x12a(i,:)) kron([0;1],x12b(i,:))];
    Z((i-1)*2+1:(i-1)*2+2,1:2*l)=[kron([1;0],z12a(i,:)) kron([0;1],z12b(i,:))];
end


[beta, T, T_J, F, part, oir] = SURPartiesRobustLagvar_IVGMM(Y, X, Z, mean_a, mean_b, RAND1, n, k, l);


% Put into final matrices
BETA(1:k, casenum)=beta(1:k);
BETA(kmax+1:kmax+k, casenum)=beta(k+1:2*k);
T_WALD(1:k, casenum)=T(1:k);
T_WALD(kmax+1:kmax+k, casenum)=T(k+1:2*k);
TJ_WALD(1:k, casenum)=T_J(1:k);
F_WALD(:, casenum)=F;
PARTIAL(1:k-1, casenum)=part(1:k-1);
PARTIAL(kmax:kmax+(k-2), casenum)=part(k:2*k-2);
DOF(:,casenum)=[k, l];
OIR(:, casenum)=oir;
clear Y X Z P beta  W F T T_J alpha ETA part mean mean_l mean_h oir;






%%%%%%%%%%%%%%%%%%%
% Compute p-values%
%%%%%%%%%%%%%%%%%%%

PVAL_TWALD=1-chi2cdf(T_WALD.^2,1);
PVAL_TJWALD=1-chi2cdf(TJ_WALD,3);
FDOF_ALL=kron([2; 1; 1],DOF(1,:));
PVAL_FWALD=1-chi2cdf(F_WALD,FDOF_ALL);
PVAL_OIR=1-chi2cdf(OIR,DOF(2,:)-DOF(1,:));





save ParticipPartiesAdjMaxLagvar_IVGMMResults BETA PARTIAL PVAL_TWALD PVAL_TJWALD PVAL_FWALD T_WALD TJ_WALD F_WALD PVAL_OIR OIR;


